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ABSTRACT 

The  problem  of  developing  orthogonal  polynomials  for  finding 
a  least-squares  best  approximation  to  a  continuous  function  over 
a  given  interval  is  solved  by  a  method  of  finding  the  mean 
derivative  of  the  function  over  that  interval.  A  parallel  method, 
using  mean  differences,  is  developed  for  fitting  discrete  data. 
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FOREWORD 


This  work  was  carried  out  as  part  of  an  independent  research 
project  entitled  Applied  Mathematics  and  Numerical  Analysis  as 
authorized  by  WEPTASK  No.  000  000  000/107 -1/R0110 101 .  The  primary 
reason  for  the  study  was  the  need  for  improved  methods  in  the 
functionalization  of  parameters  used  in  connection  with  trajec¬ 
tories  and  aiming  data.  Credit  is  due  many  of  the  present  and 
past  Naval  Weapons  Laboratory  employees.  The  following  contribu¬ 
tions  to  the  solution  of  the  problem  deserve  special  credit: 

Mr.  J.  E  T’r.rker  initiated  use  of  orthogonal  polynomials  on 
the  first  NWL  digital  computer;  Mr.  John  H.  Walker,  Mr.  Ira  V.  Wes 
and  Mr.  Larry  E.  Spangler  developed  the  technique  reported  in 
Appendix  A;  Mr.  Bruce  Bosley  and  Dr.  Hartmut  Huber  prepared  the 
tables  in  Appendix  B.  Mr.  Bruce  Bransom  made  a  study  of  the 
method  while  a  student  at  the  University  of  South  Dakota  and 
Mrs.  Jane  H.  Martina  aided  in  checking  and  editing  the  report. 
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INTRODUCTION 


This  report  is  an  outgrowth  of  the  work  reported  in  NWL 
Report  No.  1822.  Although  much  was  done  in  this  field  long 
before  NWL  Report  No.  1822  appeared,  and  outstanding  contribu¬ 
tions  have  been  made  by  others  during  the  preparation  of  this 
report,  it  is  felt  that  the  following  will  be  of  interest  to 
those  working  in  the  field: 

(1)  The  use  of  mean  derivatives  and  mean  differences  to 
tie  the  discrete  case  to  the  continuous  case  and  to 
develop  corresponding  weights  for  the  discrete  case. 

(2)  The  simple  method  of  fitting  with  certain  type 
constraints . 

(3)  The  use  of  difference  matrices  in  the  fitting  of  a 
function  of  two  independent  variables. 

Appendix  B  gives  matrices  built  from  m  values  for  making 
Jacobi-type  fits  to  equally  spaced  discrete  data.  For  higher 
degree  and  more  points  it  is  expected  that  a  high-speed  computer 
would  compute  its  own  m  values.  A  supplementary  report  is 
planned  for  giving  the  programs  to  be  used  in  computer  evalua¬ 
tion  of  w  and  m  values,  as  well  as  a  more  extensive  tabulation 
of  the  w's  and  m's  which  would  be  convenient  for  use  with  the 
desk  calculator. 

MEAN  VALUES  OF  DERIVATIVES  OVER  A  CONTINUOUS  SET 
Suppose  that  the  approximation 

f  =  Aopo  +  Aipi  +  A2P2  +  •*'  ^n  5 
where  I’k  =  Pk  (x)  .  a  kth  degree  polynomial,  a  £  x  £  b,  is 

desired.  If  equation  (1)  is  differentiated,  the  following  set 
of  equations  is  obtained: 

1  3  Ve  +  AlPl  +  a2P2  f  ••• 

i  --  A,i’J  *-  A  j  l’ *■  ... 

i"  ~  AjI’j  +  ... 


.  v>  : 


+  A  P 
n  n 


+  Vn 

*  Anpn 


(2) 


<p> 


A  p„ 

n  n 
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The  matrix  of  the  coefficients  of  the  P's  is  triangular  and  the 
system  of  equations  (2)  can  be  used  to  solve  for  the  A's  in 
terms  of  f  and  P  values . 

If  x^,  a  £  £  b,  is  substituted  in  equations  (2)  and 

the  system  solved,  a  truncated  Taylor  expansion  about  x  =  x^ 
is  obtained.  The  Taylor  expansion  is,  of  course,  in  terms  of 
derivatives  evaluated  at  a  point. 

Let  it  be  supposed  that  instead  of  using  the  k^1  derivative 
at  x  *  x^,  a  mean  of  the  derivative  over  the  set  a  s  x  s  b 
is  to  be  used.  If  w  is  an  integrAble  function,  not  changing 
signs  over  the  set  a  £  x  £  b.  normalized  so  that 

b 

J  wdx  =  1, 
a 

the  mean  kfc^  derivative  of  f  is  (omitting  the  limits  of  integra¬ 
tion,  a  and  b) 

-  /wf(k)dx  .  (3) 

If  additional  restrictions 

w(a)  ■  w'(a)  *  ...  w^k_1^(a)  *  0 

„  1X  ,  (■'♦> 
u(b)  *  w‘(b)  •  ...  (b)  ■  0 

are  made,  then  equation  (3)  can  be  changed,  Integrating  by 
partu,  to  a  form  which  might  be  easier  for  evaluation: 
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/wf^dx  *  -  Jw'f^^dx 
-  /w"f(k"2>dx 

»  (-l)k/w(k)fdx 
s  fOO 

If  w^  is  written  as  w^  ■  PQ^,  with  p  2s  0  for 

a  <  x  ^  b, 

Jpdx  *  L,  a  finite  number, 

*  a  degree  polynomial  in  x,  the  Q  polynomials 
are  orthogonal,  with  weight  p,  with  respect  to  integration  over 
the  set  a  £  x  s  b. 

This  is  shown  by 

J pQiQjdx  -  fPQ^dx  (i  <  j)  (6) 

«  Q^J)  -  0 

and 

Jpq^dx  -  Jpqj  dx  >  0.  (7) 

Special  Cases 
I.  Legendre  Polynomials 

Let  a  -  -l  and  b  -  1,  c  ■-*  a  constant  factor,  and  taka 
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cw  ■  (1  -  x2)k  , 

cw’  -  k(-2x)  (l-x2)k_1  , 

cw"  ■  k(-2)  (l-x2)k_1  +  k(k-l)  (-2x)2  (l-x2)k“2  , 

• 

and  ■  Q^,  a  pol/nomial  of  kfc^  degree. 

With  the  proper  choice  of  c,  is  the  k*-*1  Legendre  polynomial. 
This  is  the  case  for  which  p  is  constant  and  Legendre  polynomials 
are  integration-wise  orthogonal  ovez  the  continuous  set, 

-1  s  x  <1,  with  respect  to  a  constant  weighting  function. 

II.  Chebyshev  Polynomials 

Let  a  *  -1  and  b  *  1,  and  take 
cw  -  U-x2)k“1/2  , 
cw'  -  (k-1/2)  (-2x>  (l-x2)k“1-1/2  , 
cw"  -  (k-1/2)  (k-1-1/2)  (-2x)2  (l-x2)k"2_1/2  , 

+  (k-1/2)  (-2)  (1-x2)  (l-x2)k"2_1/2  , 

-  Q2(x)  (l-x2)k’2-,/2  , 

and  v0*)  m  q^d-x2) _1^2  . 

With  the  proper  choice  of  c,  Qk  is  the  k**1  Chebyshev  polynomial. 
The  weighting  function  ir  •  constant  times  (l-x2)"^2. 


<> 


III.  Jacobi  Polynomials 
For  a  ■  -1,  b  ■  1: 

cw  -  (l-x)0'**  <H*>0*  ,  a  >  -1,  0  >  -1, 
cw'  -  (O  +k)  (-1)  (1+x)  (l-x)®4*’1  (l+x)^*"1  , 

+  (/3+k)  (1)  (1-x)  (l-x/S+k"1  (l+x)^+k'1  . 

-  Q1(l-x)fl*+k”1  (l-hc)^"1  . 

cw"  -  qj(l-x)  (1+x)  (l-x)**^"2  (l+x)^*'2  , 

+  Q1(«+k.l)(-l)(l+x)(l-x)a+k‘2  (l-x^k'2  , 

+  <^(0  +k-l)  (1)  (1-x)  (1-x)**-2  (l+x>fl+k’2  , 

-  (l-x/*+k"2  (l-hc)^k"2  , 

4. 

and  -00  •  Qk(l>x)”  (1+x)^  . 

With  the  proper  choice  of  e ,  is  the  kth  Jacobi  polynomial.. 
The  weighting  function  ia  a  constant  times  (1-x)0  (i+x)^  . 

Note  that  for  a  »  0  *  0,  the  Jacobi  polynomials  become  the 
Legend  t".  polynomials;  for  a  •  £  »  -1/2,  they  bee  or*  the 
Chebvihev  polynomials . 


The  Legendre,  Chebyshev,  and  Jacobi  Polynomials  are  orthogonal 


over  the  interval  -1  £  x  £  1.  To  find  polynomials  orthogonal 
over  the  arbitrary  finite  interval,  a  £  x  £  b,  a  transformation 
in  the  independent  variable  can  be  made  to  bring  the  problem  to  the 
forms  already  studied;  or,  the  Jacobi  form  for  w  can  be  taken  as 

cw  “  (x-a)  ^  (b-x)C*+k 

IV.  Hermite  Polynomials 

For  the  infinite  interval  -oo  <  x  <  oo,  with  c  a  constant 
factor,  one  can  take 

2  2 

cw  ■  exp(- ax) 

cw'  **  -2  a2x  exp(-a2x2) 

cw,!  *  (-2a2  +  (-2a 2x)2)  exp(-a2x2) 


and  w^  =  Q,t  exp(-a'‘x2)  . 

is,  with  the  proper  choice  of  c  and  a,  the  ktn  Hermite 
polynomial,  th«  weighting  function  being  a  constant  times 
exp(-  a2x2) . 

V.  Laguerre  Polynomials 

For  the  semi -infinite  interval  Os  x  <  oo  one  can  take 


b 


=  xk  exp(-dx) 

=  (k»ax)xk^  exp(-ax) 

=  [-ax  +  (k-ax)  (k-i-ax)]xk_1  exp(-ax) 
=  Q2xk“2  exp  (-ax) 


and  w<k>  =  exp(-ax)  . 

With  the  proper  choice  of  c  and  a  ,  is  the  Laguerre 
polynomial;  the  weighting  function  is  a  constant  times  exp '-ax). 

The  literature  gives  a  more  complete  discussion  of  orthogonal 
polynomials  and  orthogonal  functions  (see  reference  (1)).  It 
should  be  pointed  out,  that  while  use  of  orthogonal  polynomials 
for  the  P's  in  equation  (1)  leads  to  a  diagonal  matrix  of  coeffic¬ 
ients  of  the  A's  in  the  set  of  linear  equations  obtained  by  taking 
mean  ktl1  derivatives,  any  desired  form  for  the  P's  can  be  used  and 
the  resulting  matrix  of  coefficients  of  the  A's  will  be  no  worse 
than  triangular.  Further,  any  general  non-negative  w  function  can 
De  used  in  equation  (3)  for  finding  the  mean  kc^  derivative,  and 
the  resulting  matrix  of  the  coefficients  of  the  A's  will  be  no 
worse  than  triangular. 

LEAST  SQUARES  FITS  IN  TERMS  OF  JACOBI  POLYNOMIALS 


Continuous  Data 


It  will  be  assumed  that  the  fit  is  over  a  f  to  interval 
a  ■  x  b,  and  that  a  linear  transformation  has  been  made  in  the 


argument,  u  =  L(x)  ,  in  such  a  way  that  with  the  function  expressed 


as  f(u)  the  fit  is  made  over  the  interval  -1  s  u  s  1, 
The  orthogonal  polynomials  have  the  property  that 


/  pCujQ^uX^CuJdu  =  0,  i  \  j, 


with  p(u)  >eing  the  weighting  function. 

The  weighted  least-squares  fit 

f  (u)  =  AqQ0  (u)  +  AjQ^  (u)  +  A2Q2  (u)  +  ...  +  A^Qj. 

is  made  by  defining 

e(u)  =  f(u)  -  AqQ0(u)  -  AjQ^u)  -  A2Q2(u)  -  ... 


and  minimizing 


E  ■  j  p(u)e2(u)du  . 

-1 


Setting 


«  ° 

caT 


leads  to 


(8) 


(u)  (9) 


'  \V“> 
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Discrete  Data 


Suppose  f(x)  I3  given  for  a  discrete  set 
a  =  Xq  <  x^  <  X2  ...  <  xm  =  b ,  and  that 


of  values  of  x, 
the  transformation 


2x;-a-b 

u  =  “b^T” 

has  been  made  so  what  the  corresponding  arguments,  in  u,  are 


(12) 


u0  =  'l  <  ui  <  u2  <  <  um  =  l- 

Replace  the  integrals  in  equations  (10)  by 

1 

/  P(u)Qk(u)f  (u)du 

(ui  -  [  Wf<V  +  VVf(ui>] 

+  (u2  •  ul>  +  Qk<u2)£(u2>] 


(“» •  Vl>  (!'■  V,)  [VVll«Vl>  +  Wf<u«>]  '  <»> 

N  m  m-1'  u  J 
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1 

/  PQ<  (u)Qk(u)du  ^ 

1  J 

(,,i  •  v  <^^[vvVuo> +  V^w' 


+  • 


+  (u  -  u  ,) 
'  m  m-1 


(UI(.'u’m))'  k<".-i>\<vi> +  V"*)Q«(v>]  (14> 

v  m  m-1'  L  j  j 


where 

us 

I (r  ,s)  =  i  /  P(u)du  .  (15) 

2  ur 

Note  that  this  is  the  result  of  replacing  p(u)  in  each  interval 
u  *  u  us  by  its  mean  in  that  interval  and  then  integrating 
by  the  trapezoidal  rule. 

In  matrix  language,  let  the  set  of  equations  obtained  by 
substituting  i;  3  u^,  u^,  Uj,  ....  un  in  equation  (9)  be 


Q  A  =  F 


(lb) 


11 


where 


12 


superscript  *  indicatirg  transpose  and  D  being  a  diagonal 
"weighting"  matrix 


.1(0,1)  0 


•  •  •  0 


1(0,2)  0 


1(1,3)  ...  0 


1  (m-1  ,m) 


Special  Cases 
A,  Legendre  Polynomials 


In  this  case,  P( u)  =  1  and  I(r,s)  «  A  (u  -  u  ) 

2  s  r' 

B.  Chebyshev  Polynomials 


In  this  case,  P(u)  =  (1  -  t,2) 


2s  2 


Ur.s)  =  i  [ arcsine  ug  -  arcsine  j 

XOte  U,at  thc  matrix'  n  Q.  in  equation  (17)  will  not  be 
diagonal,  but  will  be  near  enough  diagonal  in  form  to  give  a  good 
inverse . 


1) 


If  a  P-matrix  is  built  in  the  same  manner  as  the  Q-matrix 
was ,  the  least  squares  solution  of  P  ^  ^  can  be  found  by 

solving 

Q*  D  P  =  Q*  D  F  .  (18) 

The  proof  of  this  can  be  constructed  from  the  fact  that  every 
PK(u)  is  a  linear  combination  of  orthogonal  polynomials  of 
degree  K,  K-l,  —  1,  0.  The  Q*  D  P  matrix  will  be  almost 
triangular. 

Needed  values  of  Legendre  polynomials  can  be  computed  using 
P0(u)  =  1 

P^u)  =  u 

P2(u)  =  i  (3u2  -  1) 

Pr+i<u)  =  [(2r+l)uPr  (u)  -  rPr_1(u)J 

Needed  values  of  Chebyshev  polynomials  can  be  generated  using 
T0(u)  *  1 

Tj (u)  *  u 

T2(u)  *  ?u2  -  1 

Tr+1(0  -  2uTr(u)  -  Tj-.^u)  . 
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LEAST  SQUARES  FITS  TO  DISCRETE  DATA  OBTAINED 
BY  MEAN  DIFFERENCES 

Let  the  form  desired  for  the  fit  to  (approximation  of)  a 
function  be 


If  f(u)  is  known  for  a  discrete  set  of  values  of  u  (greater 
than  (n+1)  in  number)  and  substitution  of  these  u  values  are  made 
in  equation  (!')  the  result  can  be  written 


with  each  row  of  the  U  matrix  being  the  row  matrix  in  equation  (1') 
evaluated  for  one  of  the  u  values  for  which  f(u)  is  known  and  the 

~N  — * 

vectors  A  and  F  being  defined  as  in  equation  (17) . 

Suppose  chat  multipliers  can  be  found  such  that  if  these 
multipliers  be  used  as  components  of  a  vector,  the  dot  product  of 
this  vector  and  the  F  vector  is  a  constant  times  the  mean  Uth 
difference  ot  f.  A  matrix  M  can  then  be  built  with  the  first  row 
be  inn  the  multipliers  tor  t  fading  the  mean  zeroth  difference, 
the  second  row  being  the  multipliers  for  finding  the  mean  first 

difference . and  the  (n+lht  row  being  the  multipliers  for 

finding  the  mean  nth  difference.  Then  premu  it  i  plv  tng  !<v  M,  toe 

1  ■ 


the  equation 


M  U  A  -  M  F 

is  obtained. 

But  (see  reference  (2)) 

M  U  -  T 

with  T  a  nonsingular  triangular  matrix.  Therefore, 

A  -  (MU)’1  M  F  (19) 

-1 

«  T  1  M  F 

It  should  be  noted  that  construction  of  the  multipliers  for 
unequal  spacing  is  much  less  attractive  than  is  the  case  for  equal 
spacing,  and,  in  general,  is  to  be  avoided  unless  the  particular 
grid  of  the  unequal  spacing  is  to  be  used  repeatedly. 

Identification  of  fits  thus  obtained  with  least  squares  fits 
can  be  made  by  following  a  path  parallel  to  that  used  in  developing 
the  Jacobi  polynomials  by  finding  mean  derivatives. 

Let  p  1  be  a  positive  definite,  synsnetric,  weighting  matrix 
and  Q  be  a  matrix  whose  (k+l)st  column  is  built  from  evaluations, 
at  a  discrete  set  of  arguments,  of  a  k^  degree  polynomial,  Qk , 
the  %'s  being  constructed  in  such  a  way  that 

Q*  H~l  -  M. 

<!>  ■  fl ' 1  Q  -  H  Q  -  0 , 


In 


Then 


i 


a  diagonal  matrix,  and  the  Q  polynomials  are  suimnattonwise 
orthogonal  with  respect  to  the  given  weighting  matrix.  The 

proof  is  that  for  j  >  i  the  jth  mean  difference  of  an  ith  degree 
polynomial  is  zero. 

The  weighted  least  squares  solution  of  V  t  •  t  ia  the 
solution  of 

U*  M  ^  I)  A  *  U*  f\ 

But 

P0  *  bQC^0  ’ 

Vb0lVVl  * 

Pn  *  b0n%  +  binQi  +  •  •  •  bnnQn  * 
or  U  *  Q  B, 

with  B  «  nonsingular,  triangular  matrix.  Therefore 

u*  m'1  vt  -  u*  41-1  f 

can  be  written 

B*  Q*  4'1  U  A  .  B*  Q*  fT  l  $ 

or  B*  H  l)  A  »  B*  M  F 

and  since  B  is  nonsingular 

M  U  A  *  H  F  % 


MEAN  VALUES  L?  DIFFERENCES  OVER  A  DISCRETE  SET 
The  problem  of  finding  the  mean  ktn  difference  from  a  set  of 
functional  values,  defined  over  a  discrete  set,  has  been  discussed 
in  reference  (2).  A  development,  similar  to  that  indicated  in 
equation  (4) ,  is  imnediately  seen  by  writing  for  r+1  equally 
spaced  data,  with  k  =  2  (see  reference  (2)),  the  weighted  mean 
second  difference  equals 

(/j2y)  =  [Ml(y0  ‘  2yl  +  y2)  +  W2(yl  '  2y2  +  y3> 

+ . -  - 

+  wr-l(yr-2  -  2yr-l  +  yr>  ]  * 

In  matrix  language  this  becomes 


ON 


L. 


3 


19a 


>0 

yl 


Reference  (3)  makes  the  steps  to  equation  (20)  parallel  to  those  of 
equation  (4)  by  using  summation  by  parts.  For  finding  the  mean 
difference,  the  steps  leading  to  equation  (20)  show  that  the 
m’s  are  (-1)^  times  the  kfch  differences  of  the  set 

[000  .  .  .  0wjW2  .  .  .  •  •  •  00] 

having  k  zeros  at  each  end  (reference  (2)). 

The  Gram  polynomials,  which  are  the  m's  for  equally  spaced, 
equally  weighted  data,  have  been  discussed  in  many  sources.  A 
partial  tabulation  may  be  found  in  reference  (2)  and  in  many 
other  published  works  on  least  squares  fitting.  An  interesting 
game,  for  any  computer,  is  to  see  how  iur  one  can  go  in  number 
of  equally  spaced  points  and  degree  of  polynomial  with  the  m's 
represented  exactly  by  the  number  of  digits  used  in  single 
precision  computation.  (See  Appendi>  A.) 
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When  this  limit  is  exceeded,  the  w's 


can  be  cut  in  number  of  digits,  and  these  approximate  w's  used 
to  find  the  m's  within  the  capacity  of  single  precision  arithmetic. 

By  considering  the  weighting  function  for  getting  the  mean 
k^1  derivative  which  leads  to  the  development  of  the  Jacobi 
polynomials,  it  is  possible  to  construct  a  corresponding  set  of 
weights  for  finding  the  mean  kth  difference.  Take 

a+k  fi+k 

cw  =  (1-x)  (1+x) 

k(l+S)  k(i+£) 

-  (1-x)  (1+x) 

k  (1+  ®  )  (1+  £  ) 

“  [(1-x)  ]  [(1+x)  ]  ,a>-l,/ff>-l. 

When  this  form  is  compared  with  the  a  ■  0,  f}  *  0  case,  which 
leads  to  the  Legendre  polynomials,  to  the  weights 


used  for  equally  spaced,  equally  weighted,  data  in  finding  the  mean 
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k^  difference  corresponding  to  the  Gram  polynomial,  it  is 


seen  that  the  set  of  weights. 


"/ r+l-i  \  ~]  *  r/  k-l+i  \  "I 

i-L(  *  )J  L(  „  )J  •k>0  (22> 


(1+2) 

k 


(1+  f  ) 

k 


can  be  used  for  equally  spaced  discrete  data  in  finding  fits  corres¬ 
ponding,  roughly,  to  the  least-squares  fits  made  with  Jacobi  poly¬ 
nomials  over  a  continuous  set.  It  should  be  noted  that  for  k  *  0 
the  w's  become  indeterminate;  for  this  case,  rounded  values  of 
the  elements  of  the  diagonal  matrix  D,  in  equation  (17),  can  be 
used. 

In  particular,  corresponding  to  the  Chebyshev  case,  if 

a  =  =  -  1/2, 


a  -  ) 

a  -  JL 

rr1-1) 

1  Zk  r/ k-i+i  > 

il  2k 

L\  k  / 

J  L\  k  ) 

(2ti) 

1  r+l-i  \  "I  2k 

ri+i) 

.  k 

L\  k  / 

\  k  /J 

These  w'c  can  easily  be  constructed  from  those  tabulated  in 
reference  (2).  With  this  in  mind,  equation  (23)  can  be  written 


(cwj). 


/ 2k - 1  \ 


k  >  0 
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The  T  subscript  indicates  fitting  with  Chebyshev  type  weights  and 
the  G  subscript  indicates  fitting  with  Gram  (or  Legendre)  type 
weights.  The  (cw^)^,  values  of  equation  (24)  are,  in  general, 
irrational.  Therefore,  they  must  be  rounded  before  differences 
are  taken  to  find  the  exact  m's  corresponding  to  the  approximate 
w's.  In  general,  if  three  significant  digits  are  retained  in 
the  first  (smallest)  nonzero  w,  a  satisfactory  set  of  m's  can 
be  found. 

The  indeterminate  form  for  equation  (22)  when  k  =  0  and 

?  9 

a  +  p  >  0  can  be  avoided  by  writing,  for  the  continuous 
case , 

k-1  1+a  k-l  1+0 

cw  -  (1-x)  (1-x)  (1+x)  (1+x)  ,  k  >  0, 


cw  = 


(1-X)a  (l+x)0  , 


0. 


Equation  (22)  is  then  replaced  by 

1+0 


/r+l-i\  / r+2-k-i  \  /k-l+i \  /i\ 

“i-U-iJl  i  )  U-iJlJ 


l+/> 


and 


,for  k  >  C , 
(22a) 


ewi  =  (r  +  2  -  k  -  i)01  (i)^  ,  for  k  «*  0 

1+a  1+  0 

The  (  ,  )  (  i  )  ,  (r  +  2  -  k  -  if  ,  (i)^  arc  to 

be  evaluated  only  when  r=2-k-i>0  and  i  >0. 

Note  that  for  all  cases  corresponding  to  Jacobi  type  fits, 
1  *  d  >  0  and  1  +  0  >0. 
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Reference  (4)  discusses  power  series  approximations  of 
functions  having  the  d'Alembert  characteristic.  Special  cases 
are: 

3  5 

fl  =  AjX  +  A-jx  +  A^x  +  . . .  , 
f^  =  A^x4  +  A&x6  +  Agx8  +  . . . , 
f  5  =  AgX^  +  k-jx'  +  A^x^  +  . . . 

These  are  functions  having  a  set  of  derivatives  specified  at  a 
point.  The  m's,  computed  from  a  suitable  set  of  w's,  can  be  used 
for  making  polynomial  approximations  of  functions  of  this  type. 

Suppose  a  transformation  has  been  made  in  the  independent 
variable  so  that  u  ~  0  for  the  value  of  x  for  which  certain 
derivatives  are  specified.  The  complete  function  can  be  approxi¬ 
mated  by 

f(u)  =  Aq  +  A^u  +  A2-u"  +  .  .  .  +  A^u  .  (25) 

Suppose  f(0),  f'(0),  f " '  (0) .  .  .  .  f (0)  are  specified.  From 


f  (0)  =  An 
f'(0)  =  \{ 


it  is  seen  that  the  problem  is  that  of  fitting 
g(u)  *  f(u)  -  Aq  -  A1u  -  ...  -  Arur 

*  A2 u^  +  A^u^  +  .  .  .  A^u"  ,  (26) 

and  that  the  unspecified  A's  are  to  be  chosen  in  a  way  consistent 
with  the  least-squares  weighting. 


Substituting  in  equation  (26) , 


If  the  matrix  M  is  constructed  by  making  the  fJ_it  row  the 
m's  for  k  *  2,  the  second  row  the  m's  for  k  »  4 ,  ....  the  last  row 
the  m's  for  k  “  n,  premultiplying  both  sides  of  equation  (27)  by 
M  will  give  the  required  number  of  equations  for  solving  for  A2 , 

A/, ;  ....  An,  .  nd  the  matrix  of  the  coefficients  of  the  A's  will  be 
a  triangular,  nonsingular .  matrix.  The  method  can  be  used  when  any 
of  the  A's  of  equation  (1)  arc  specified. 


FITTING  A  FUNCTION  OF  TWO  INDEPENDENT  VARIABLES 


Let  it  be  supposed  that  values  of  the  function  are  known  for 
rectangular  grid,  equally  spaced  in  both  independent  variables. 
Make  linear  transformations  in  the  independent  variables  so  that 
the  new  independent  variables  become  u  and  v. 

Approximate  f  by 

fPn(v)l 


’I/" 


A  *  A 


-0(u) 

P^u) 

p2  (u) 

...  P^(u) 

^0 

A01 

*02  ' 

A10 

A11 

A12 

A20 

• 

A21 

A22 

•  i 

k2V 

• 

^/O 

V, 

•  < 

A//2 

k/uv 

P^v) 
p2  (v) 


pi/(v) 


Cons  truct 


(MT)‘l  MVI  -  (T*1  M) 


tor  the  number  ot  u  v.ilucs  in  the  grid  and  the  decree  in  u  desired, 


as  In  equation  (19) .  In  the  like  manner  construct 


so  that  eight  digit,  ten  digit,  fifteen  digit,...  approximations 
will  result  with  corresponding  round-off  errors  in  the  A  matrix. 

END  POINT  SMOOTHING  AND  PREDICTION 
FROM  EQUALLY  SPACED  DISCRETE  DATA 

The  w  function,  cw  ■  x^e"0*  ,  for  finding  the  mean 
derivative  over  the  semi-infinite  interval,  0  i  x  £  oe,  which 
leads  to  the  Laguerre  polynomial,  indicates  a  w  form  that 
can  be  used  to  solve  the  problem  of  end-point  smoothing  and 
prediction  from  equally  spaced  discrete  data. 

Suppose  that  the  function,  y,  is  known  for  u  ■  0,  1,  2,  3,.. 

Take 


Although  the  Laguerre  polynomials  are  based  on  the  se^i- 
infinite  interval,  0  £  x  <  oe,  and  the  above  indicated  u-set 
would  go  to  infinity,  in  actual  practice  a  finite  number  of 
points  would  be  used. 

For  k  >  0,  the  w's  would  be  modified  by  replacing  the  last 
few  w's  with  multiples  of  those  used  in  Jacobi-type  fittings. 
This  "splicing"  would  done  so  as  to  bring  a  gradual  descent 
of  the  w's  to  zero  ami  avoid  unwanted  magnitudes  of  the  result¬ 
ing  m's. 
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A  set  of  w's  suitable  for  the  solution  of  this  problem  can 
also  be  based  cn  the  continuous  data  w's  for  Jacobi  polynomials. 
ine  set  is  formed  by  choosing  fi  £  0  and  Ct  >  0  in  equation  (22) 
for  left-end  smoothing;  for  right-end  smoothing  reverse  the  order 
of  the  nonzero  w's. 


One  application  of  end-point  smoothing  is  the  problem  of 

finding  the  time  of  motor  separation  from  observed  positions  of 

a  rocket-propelled  missile.  The  difference  between  the  right- 

hand  smoothed  estimate  of  (acceleration)^  and  the  left  hand 

smoothed  estimate  of  (acceleration)^^  would  be  expected  to 

have  maximum  magnitude  when  separation  takes  place  between  t^ 

and  d .  , . 
l+l 

Another  application  is  that  of  aiming  at  a  moving  target. 

In  this  case  the  best  possible  estimate  of  present  position  vector 

and  its  derivatives  is  desired. 

MID-POINT  SMOOTHING  FOR  AN  ODD  NUMBER 
OF  EQUALLY  SPACED  DATA  POINTS 


The  w  function,  cw  *  e 


,  for  finding  the  mean  kc 


derivative  over  the  infinite  interval  -Ct  <  x  <  oc  can  be  used 

for  this  problem.  For  discrete  data  the  first  estimate  of  the  w's 

1  2 

would  be  generated  from  cw  -1  e*0^"  with  u  *  -m,...-3,  -2,  -1, 


i).  1  ,  2 


i . m  and  linal  v’  •  can  be  obtained  hv  replacing  at 


either  end  with  w's  from  a  Jacobi-type  fit.  The  need  for  this 
modification  is  not  so  severe  as  in  the  case  of  end-point 
smoothing. 

This  problem  can  also  be  solved  by  using  equation  (22)  and 
setting  Ct  =  0  >  0. 
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APPENDIX  A 


COMPUTATION  OF  w  AND  m  VALUES  BY  HIGH  SPEED  COMPUTER 


For  computational  purposes,  every  number  whose  factorial  must  be 
evaluated  is  represented  in  its  prime  factor  form.  This  is  necessary 
in  order  to  conveniently  represent  the  w's  in  the  equation, 

(n  -  i)  1 


(k  -  i  +  i) : 
c(^i  ‘  k!  a  -  d: 


k!  (n  -  i  -  k) ! 


(1) 


where 

n  =  number  of  data  points 

and 

k  =  the  order  of  difference. 

If  the  computations  were  performed  in  decimal  form,  the  resulting 
value  of  (wk)i  might  well  be  expected  to  exceed  the  IBM  7030's  maximum 
word  length  of  14  digits.  The  product  of  the  numerators  of  equation  (1) 
may  be  formed  by  adding  the  exponents  of  the  corresponding  prime  factors 
of  each.  Likewise,  the  product  of  the  denominators  may  be  formed  in  the 
same  way.  Division  may  be  accomplished  by  subtraction  of  the  exponents 
of  the  corresponding  prime  factors  of  each,  resulting  in  a  value  for 
cCwj^  which  is  the  product  of  15*  prime  factors  raised  to  various 
powers . 


*  The  number  15  was  chosen  because  in  the  computer  program  15  prime 
factors  are  needed  to  represent  the  factorials  of  all  the  numbers 
from  1  to  50.  They  are:  2,  3,  5,  7,  11,  13,  17,  19,  23,  29,  31, 
37,  41,  43,  and  47. 


1 


After  the  (n  -  k)  values  of  c(w^)^  are  computed,  the  (wk)/s  are 
reduced  to  their  lowest  terms.  Then,  evaluation  of  the  decimal  equiva¬ 
lent  of  each  prime  factor  (w^)^  can  be  made.  This  is  accomplished  by 
first  obtaining  5  decimal  numbers,  each  number  computed  by  forming  the 
product  of  3  of  the  15  quantities  in  the  prime  factor  representation 


of  (wk)^.  The  following  illustration  may  clarify  the  procedure: 
.  „  Q1  ^3  _Q4  .„Q6  .  Q7  .  Q8 


(wk)i  =  2  -3 


11  ’  13  •  17  •  19  •  23 


^10  ^11  ^12  ^13  ^14  Q15 

29  *  31  *  37  •  41  •  43  •  47 


where  Qm  is  an  integer  or  zero  (m  =  1,  2  .  .  .  15) • 


Ql  Q6  Qh  Q2  Q7  Q12 
(wR)i  =  (2  -13  *31  )  *  (3  *17  -37  )  .  .  .  . 

Q5  Q10  Q15 

(11  *  29  1U  •  47  ) 


Rj.  ’  R2  •  •  *  •  Rft  >  Cm  -  1 ,  2 ,  . . .  5) , 

where  each  expression  in  parentheses  (and  therefore.  .  .  .  R^) 

is  a  decimal  number. 

At  this  point,  each  R^  is  examined  and  in  every  case  the  winner  of 
places  to  the  left  of  the  decimal  of  each  is  stored.  These  "magnitude 
tags"  are  then  summed  and  the  result  is  an  indication  of  the  magnitude 
01  each  (wk) j .  If  the  sum  of  the  "magnitude  tags"  is  less  than  15, 
multiplication  of  the  R^j ' s  is  performed  and  a  («k) ^  is  obtained.  If  the 
sum  oi  the  "magnitude  tags"  is  13  or  more,  R^  is  repeatedly  divided 


by  10  until  the  sum  of  the  "magnitude  tags"  indicates  that  the  new 
(reduced)  (w^) ^  can  be  represented  in  less  than  15  digits.  The  R^'s 
are  multiplied  and  each  (w^)^  value  obtained  is  truncated  to  an 
integer. 

Thus,  having  computed  the  the  next  step  is  to  obtain  the 

m's.  This  is  done  by  evaluating  the  matrix  equation,  ra  ■  [b]  , 


and 


3 


with  k  zeros  at  each  end. 

That  each  (j\)  t  may  be  contained  in  less  than  15  digits  is 
assured  by  employing  another  "magnitude  tag"  before  each  multipli¬ 


cation  . 


APPENDIX  B 


TABLES  OF  TRANSPOSE  OF  M  ANT)  INVERSE  OF  MU 


These  tables  are  designed  for  finding  the  weighted  least  squares 
fits  in  the  form 

2  i  4  5 
y  ■  A_  +  A,u  +  A„u  +  A,u  +  A,  u  +  A,u 
'  0  1  2  3  4  5 


to  2s  +  1  data  points,  u  taking  the  values 

-s,  (-s  +  1),  ...  -2,  -1,  0,  1,  2,  ...  s  -  1,  s. 

The  w's  used  in  constructing  the  M  matrices  were  constructed  by 
ordering  the  distinct  non-zero  computed  w's  so  that 

Wj  £  v?2  <>  w3  <;  •  •  .  wm 

and  choosing  the  q  from  the  integers  1,  2,  3,  ...  N,  N  s  25, 
which  makes 


a  minimum,  being  the  integer  nearest 

The  computed  w. ,  w,,  ...  w  values  were 
i  m 


( 

wi 


then  replaced  with 


Pi  *  q’  P2  ’  P3  ’  •”  V 

,'ig'n  '  -e  c'-~>sen  so  is  to  make  the  resulting  fits  correspond, 
cilv.  to  least  squares  fits  using  Jacobi  polynomials  with  weighting 
* 1 '  at  ed  hv  the  given  (  a  ,  f}  )  values.  ihus  (O.O^  indicates  that  a 
degree  *  it  is  being  made  with  toe  "Jacobi"  a  .  ft  1  equal  to 


l 


(0,0).  This  is  the  least  squares  fit  with  uniform,  or  equal, 
weighting.  (-1/2,  -1/2)5  indicates  that  a  fifth  degree  fit  is  being 
made  with  the  "Jacobi"  (Ct,0)  equal  to  (-1/2,  -1/2).  This  corresponds 
to  the  Chebyshev,  or  minimax,  type  fit,  and  is  a  least  squares  fit 
with  heavy  weighting  at  the  ends  and  smaller  weighting  at  the  center 
of  the  interval.  (1/2,  1/2)5  indicates  that  a  fifth  degree  fit  is 
being  made  with  the  "Jacobi"  (a  ,  fi)  equal  to  (1/2,  1/2).  This  is  a 
least  squares  fit  with  heavy  weighting  at  the  center  and  smaller  weight 
ing  at  the  ends  of  the  interval. 

Since  the  u  values  and  the  weighting  chosen  have  symmetry  about 
the  mid  point,  or  center,  each  (MU)  ^  takes  the  form. 


(MU) 


-1 


"00 


0 

0 

0 

0 

0 


'11 


0 

0 

0 

0 


'02 


22 


13 


33 


"04 


24 


44 


0 

C15 

0 

C35 

0 

C55 


I.et  F  be  the  vector  defined  for  equation  (17)  with  uQ  ■  -s, 

u j  =  -s  +  l,  ...  and  -  s,  and 

_\ 

bo  the  first  column  of  , 


m^  be  the  seccnd  column  of  M* 


m^  be  the  sixth  column  of  M*  . 


The  A's 
function  arc 


in  the  least  squares  polynomial  approximation  of  the 
then  given  by 


*0 

A1 

A2 

A3 

A4 


T. 

coo 

mo 

• 

?+C02 

Ok 

.  F  + 

C04 

_k 

m4 

Ok 

Ok 

y 

C11 

ml 

a 

F  +  C13 

.  F  + 

C15 

m5 

Ok 

Ok 

m 

C22 

*2 

• 

F  +  c24 

m4 

.  F 

-V 

IX 

C33 

• 

F  +  C35 

m 

5 

.  F 

C44 

-2k 

m4 

• 

F 

= 

C55 

• 

m5 

• 

-o 

F 

F 

F 


Lower  degree  fits 
from  the  right,  in  the 


can  be  made  by  deleting  corresponding  colun. 
matrices  M*  and  (MU)"*. 
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(0,0)5 


11  Points 


853 

40898 

0 

-3569 

490776 

0 

1 

3432 

0 


(-1/2,  -1/2)5 


11  Points 


1 

166 


0 


J. 

39 


27323 

1437228 


0 


1 

162  0 


40753 

931500 


0 


0 


1 

468 


-1081 

207792 


0 


1 

1500 


-287 

46000 


0 


0 


1 

5328 


0 


0 


0 


1 

5520 


-1/2)5  11  Points 


22 

-11 

11 

-11 

11 

-3 

16 

-4 

-1 

7 

-15 

7 

14 

-2 

-1 

5 

-5 

-2 

13 

-2 

-3 

5 

2 

-2 

12 

0 

_/ 

2 

2 

-3 

12 

0 

-4 

0 

10 

0 

12 

0 

-4 

-2 

2 

3 

13 

2 

-3 

-5 

2 

2 

14 

2 

- 1 

-5 

-5 

2 

16 

4 

-1 

-7 

-15 

“  / 

22 

11 

11 

11 

11 

3 

(1'2,  1/2)5 


11  Points 


(MU) 


0 


-619 

58756 


0 


7775 

940096 


0 


0 


1 

312 


-131 

11336 


110159 

14963520 


0 


0 


1 

794 


-8881 

2439168 


0 


0  0 


-637 

460416 


0  0 


1 

6144 


0 


o  0  0 


1 

21120 


(1/2,  1/2)5  11  Points 


M*  = 


9 

-9 

10 

12 

-j.3 

10 

14 

-12 

2 

15 

-9 

-5 

16 

-5 

-11 

16 

0 

-12 

16 

5 

-11 

15 

9 

—  l 

14 

12 

2 

12 

13 

10 

9 

9 

10 

-6 

9 

-9 

-2 

-5 

15 

6 

-15 

11 

7 

-4 

-14 

6 

6 

-20 

0 

18 

0 

-6 

6 

20 

-7 

-4 

14 
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APPENDIX  C 


COMPUTATION  OF  THE  ORTHOGONAL  POLYNOMIALS 
AND  WEIGHTING  MATRIX  FOR  GIVEN  MULTIPLIERS 


In  reference  (2)  a  method  for  going  from  the  weighting  matrix 
and  the  corresponding  orthogonal  polynomials  to  the  w's  and  m's 
was  indicated.  Since  rounding  of  the  v's  occurs  in  tha  construction 
of  a  convenient  set  of  m's  it  is  the  purpose  of  this  appendix  to 
show  a  method  of  construction  of  the  orthogonal  polynomials  and  the 
weighting  matrix  from  che  m's. 

Suppose  there  are  r  +  1  discrete  points, 

u  *  the  vector  whose  components  are  the  nc**  powers  of  the 
grid  point  arguments, 

*  the  vector  whose  components  are  the  multipliers  for 
finding  the  mean  difference, 

Qy  »  the  vector  whose  components  are  evaluations  of  the 

orthogonal  polynomial  for  the  grid  point  argunents, 
M  *  »  the  (r  +  1)  x  (r  +  1)  weighting  matrix, 

Tc  construct  find  lctst  squares  (k  -  l)st  degree  fit,  using 

Unown  m's,  to 

'‘k  ”  Ko%  +  *  •••  +  S-l^k.l  + 

If  ^  is  this  (k-l)st  a  -iSev  polynomial  fit 


i 


Take  all  components  of  to  be  1,  construct  Q  ,  C^,  Qr 

and  let 


Q  = 

M*  - 


J>  A  A 

%  Ql  ^2 


t 


-2k 

mo  ml  *2 


Q  and  M*  are  (r  +  1)  x  (  r  +  1)  nonsingular  matrices.  Prom 
reference  (2) . 


M*  =  M  *Q,  and  therefore 

M*Q  ^  \  the  required  weighting  matrix. 

Note  that  the  m's  obtained  from  rounded  w's  give  a  true  least 
squares  fit  with  respect  to  a  "rounded"  weighting  matrix. 
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